8. ジョルダン標準形とスペクトル集合

ジョルダン標準形

code: jordan.py

from sympy import *

from numpy.random import seed, permutation

A = diag(1, 2, 2, 2, 2, 3, 3, 3, 3, 3)

seed(123)

for n in range(10):

P = permutation(10)

B = Lambda(S('lmd'), A - S('lmd') * eye(10))

x = Matrix(var('x0, x1, x2, x3, x4, x5, x6, x7, x8, x9'))

y = Matrix(var('y0, y1, y2, y3, y4, y5, y6, y7, y8, y9'))

z = Matrix(var('z0, z1, z2, z3, z4, z5, z6, z7, z8, z9'))

a1 = x.subs(solve(B(1) * x))

a2 = x.subs(solve(B(2) * x))

b2 = y.subs(solve(B(2) * y - a2))

a3 = x.subs(solve(B(3) * x))

b3 = y.subs(solve(B(3) * y - a3))

c3 = z.subs(solve(B(3) * z - b3))

v0 = a1.subs({x9: 1})

v1 = b2.subs({y6: 1, y7: 0, y8: 0, y9: 0})

v2 = B(2) * v1

v3 = b2.subs({y6: 0, y7: 1, y8: 0, y9: 0})

v4 = B(2) * v3

v5 = c3.subs({z5: 1, z6: 0, z7: 0, z8: 0, z9: 0})

v6 = B(3) * v5

v7 = B(3) * v6

v8 = b3.subs({y6: 1, y7: 0, y8: 0, y9: 0})

v9 = B(3) * v8

if __name__ == '__main__':

V = v0

V = V.row_join(v)

print(repr(V.inv() * A * V))

実行結果:

code: python

Matrix([

Matrix([

Matrix([

固有値1の固有空間と一般化固有空間

https://gyazo.com/2d4cc316f1f848ae6231eb22c9e6b5ab

固有値2の固有空間と一般化固有空間

https://gyazo.com/5ca613b71641256c28d6b5d6149ea006

固有値3の固有空間と一般化固有空間

https://gyazo.com/9b81fed8e9a4ab5f529c9568a2063e01

ジョルダン標準形及びジョルダン分解の問題作成プログラム

code: jordan2.py

from sympy import Matrix, diag

from numpy.random import permutation

while True:

A = X.copy()

while 0 in A:

i, j, _ = permutation(3)

if max(abs(A)) >= 10:

break

if max(abs(A)) < 10:

break

U, J = A.jordan_form()

print(f'A = {A}')

print(f'U = {U}')

print(f'U**(-1)*A*U = {J}')

B = A - C

print(f'B = {B}')

print(f'C = {C}')

実行例:

code: python

行列のスペクトル集合

code: spectrum.py

from numpy import matrix, pi, sin, cos, linspace

from numpy.random import normal

from numpy.linalg import eig, eigh

import matplotlib.pylab as plt

N = 100

B = normal(0, 1, (N, N, 2))

Real = A + A.conj()

Hermite = A + A.H

PositiveDefinite = A * A.H

PositiveComponents = abs(A)

X = PositiveComponents

r = max(abs(Lmd))

T = linspace(0, 2 * pi, 100)

plt.plot(r * cos(T), r * sin(T))

plt.scatter(Lmd.real, Lmd.imag, s=20)

plt.axis('equal'), plt.show()

複素行列

https://gyazo.com/c151fce2f3971edaaaa1f169a710da4f

実行列

https://gyazo.com/d984bda740b6307aa359ef55ed3e8de4

エルミート行列

https://gyazo.com/82e2bf899628825d957c0c2949858047

正定値行列

https://gyazo.com/a9db4e3e99578704b6fe1b36426e0157

ユニタリ行列

https://gyazo.com/840a32bc03320cfab58f9569abc915c7

正成分行列

https://gyazo.com/c7dccc0c52d047875cd9971c548e5b1f

ノルム公式

code: norm.py

from numpy import matrix

from numpy.linalg import eig, norm

from numpy.random import normal

import matplotlib.pyplot as plt

def power(ax, m):

A = matrix(normal(0, 1, (m, m)))

N = range(50)

for i in range(m):

for j in range(m):

fig, axs = plt.subplots(2, 5, figsize=(20, 5))

plt.show()

https://gyazo.com/085359f1be61afdcd25697ac04745247

ゲルファントの公式

code: gelfand.py

from numpy import matrix

from numpy.linalg import eig, norm

from numpy.random import normal

import matplotlib.pyplot as plt

def gelfand(ax, m):

A = matrix(normal(0, 1, (m, m)))

N = range(50)

fig, axs = plt.subplots(2, 5, figsize=(20, 5))

plt.show()

https://gyazo.com/48fdb811abb6b2229f0d2c31e7dd2626